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ABSTRACT 

A  three-dimensional  blunt  body  inverse  technique  is  formulated  in  terms 
of  two-component  stream  functions  in  the  spirit  of  the  Ferri,  Vaglio- 
Laurin  axisymmetric  blunt  body  method.  In  a  similar  manner  to  the 
axisymmetric  case,  transformations  are  found  that  reduce  the  elliptic 
asymmetric  shocklayerto  a  transformed  space  in  which  shock,  stream¬ 
lines  (or  the  intersections  of  stream  surfaces),  body,  and  sonic  surface 
are  a  priori  known,  so  that  numerical  analysis  procedures  are  consid¬ 
erably  simplified.  Provided  three  equations  in  the  three  velocity  com¬ 
ponents  of  the  transformed  space  are  numerically  solved,  then  inverse 
transformations  are  given  so  that  the  corresponding  velocities,  pres¬ 
sure,  and  density  at  a  point  in  the  physical  space  may  be  determined. 
The  class  of  given  shocks  considered  are  portions  of  prolate  ellipsoids 
in  the  physical  space. 
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I.  INTRODUCTION 


In  order  to  predict  the  subsonic  inviscid  flow  field  region  about  blunted  bodies 
without  axial  symmetry,  or  axisymmetric  or  asymmetric  bodies  at  an 
angle  of  attack,  a  detailed  analysis  of  an  improperly  set  elliptic  system  of 
nonlinear  equations  must  be  solved  by  a  stable  finite  difference  scheme. 

For  axisymmetric  bodies  (at  zero  and  at  small  angles  of  attack)  many 
problems  associated  with  a  straightforward  approach  to  the  above  problem 
have  been  minimized  by  transformation  to  a  new  space  upon  whose  bound¬ 
aries  the  boundary  conditions  are  known  (Refs.  1  and  2),  and  where  velocity 
distributions  at  the  axis  and  in  the  neighborhood  of  the  sonic  line  are 
a  priori  determined  and  used  as  boundary  data.  In  addition,  the  domain  of 
influence  of  the  sonic  line  is  known  a  priori,  leading  to  uniqueness  of  the 
computed  body  for  a  given  analytic  shock;  further  reasons  for  attempting 
an  extension  of  this  method  for  three-dimensional  flows  are  obvious  from 
the  discussion  presented  in  Ref.  2. 

Previous  analysis  of  regions  of  this  type  (Ref.  3)  by  series  expansion 
methods,  with  a  Clebsch  transformation,  generally  leads  to  a  great  deal 
of  computation,  indicating  the  necessity  of  an  analog  to  the  comparatively 
simple  and  theoretically  sound  method  of  Refs.  1  and  2.  In  most  cases  of 
Ref.  3,  two  or  three  terms  of  the  expansion  (involving  the  solution  of  the 
associated  ordinary  nonlinear  coupled  systems  of  equations)  were  necessary 
for  shock  and  sonic  line  shapes. 
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At  present,  numerical  methods  of  the  kind  necessary  to  fully  understand 
error  analysis  for  the  type  of  calculations  suggested  by  the  inverse  prob¬ 
lem  have  not  been  fully  developed;  ^  however,  the  methods  of  Refs.  4  and  5 
suggest  approaches  for  improperly  set  Laplace  or  Cauchy- Riemann  systems 
with  regard  to  round-off  error  analysis,  approximation  error,  necessary 
step  size,  etc.,  for  the  linear  problem. 

Provided  the  sonic  surface  is  known  (or  the  elliptic  solution  canbe  extended 
sufficiently  past  the  sonic  surface  as  is  sometimes  done),  a  number  of 
three-dimensional  characteristics  methods  can  be  used^  (Refs.  10,  11, 
and  12)  to  continue  the  solution  into  the  hyperbolic  domain. 

In  accordance  with  this  discussion,  then,  an  inverse  method  is  developed 
in  the  spirit  of  the  Ferri,  Vaglio-Laurin  method  for  three-dimensional 
flows  about  nonaxisymmetric  zero  angle  of  attack  bodies.  A  transforma- 
tion  of  the  region  R^  (Fig.  1 )  to  a  region  R  (Fig.  2)  is  effected,  and  sub¬ 
sequent  development  of  a  suitable  stream  function  pair  leads  to  a  system 
of  three  nonlinear  equations  in  the  unknown  velocity  components  in  R, 
which  are  solvable  by  the  method  of  Refs.  1  and  2,  when  extended  to  three- 
dimensional  nets.  Solutions  in  the  transformed  space  are  then  transform¬ 
able  into  the  physical  space  for  the  actual  flow  and  resulting  streamlines. 

- -  j 

1  Reflection  methods  for  elliptic  systems  have  been  used  to  obtain  symmetric 
hyperbolic  transformed  systems  for  the  flow  about  blunted  bodies  at  a  small 
angle  of  attack  (Ref.  4).  For  such  systems,  stable  numerical  techniques 
and  error  estimates  are  to  be  found  (for  example  Ref.  5  through  9)  in  con¬ 
trast  to  the  "marching"  method  for  improperly  set  elliptic  equations. 

^  Many  additional  comments  concerning  the  numerical  and  theoretical  methods 
of  Ref.  8  are  presented  in  Vol..  II  of  the  report  (in  preparation). 
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In  that  the  axis  flow  is  not  completely  a  priori  determined  from  the 
resulting  system  of  equations,  it  is  recommended  that  the  bounding 
surfaces  be  taken  as  the  shock,  sonic  surface,  and  body  approximation 

■"U 

surface  in  the  R  space  for  purposes  of  computation.  From  a  computa- 
tional  standpoint,  round-off  errors,  say  on  the  initial  surface,  propagate 
exponentially  and,  therefore,  obviate  methods  of  filtering  high  frequencies 
for  stable  solutions;  alternatively,  for  a  high  accuracy  difference  scheme, 
only  a  few  steps  are  necessary  across  the  shock  layer  so  that  growth  of 
round-off  error  is  lessened.  For  known  flow  boundaries,  the  latter 
alternative  leads  to  difference  equations  always  satisfying  the  sonic 
surface  given  boundary  condition.  From  Ref.  2,  from  the  fact  that  stream 
surfaces  are  a  priori  known,  a  lower-order  difference  scheme  is  neces¬ 
sary,  since  the  order  of  the  system  has  been  reduced.  Although  details 
of  the  numerical  analysis  will  not  be  dealt  with  in  detail  here,  it  is  possible 
to  see  that  many  of  the  numerical  simplifications  of  Refs.  1  and  2  carry 
over  to  the  method  to  be  indicated.  Although  solutions  may  be  computed 
by  a  variety  of  methods  extended  from  their  axisymmetric  counterpart, 
the  method  of  Refs.  1  and  2  appears  to  be  singular  in  the  use  of  known 
information  on  the  axis  and  sonic  surface,  leading  to  a  more  theoretically  ’ 
sound  and  numerically  simple  method  in  the  extension  to  a  higher 
dimension. 

For  a  specified  asymmetric  portion  of  the  shock  surface  given  over  the  sub¬ 
sonic  and  transonic  regions  (here  taken  as  a  portion  of  a  prolate  ellipsoid), 
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the  method  indicated  in  this  report  may  be  applied  to  determine  the 
three-dimensional  shock  layer  in  detail  and  also  the  body  shape. 

Knowledge  of  the  aerodynamic  characteristics  (pressures,  velocities, 
etc.)  in  the  nose  region  of  a  particular  "lifting"  surface  may  be  deter¬ 
mined  by  this  method,  and  inputs  to  a  downstream  characteristics- 
routine  established.  In  this  manner,  total  aerodynamic  analysis  of  fully 
three-dimensional  lifting  bodies  may  be  theoretically  performed  in  order 
to  obtain  design  criteria  for  such  bodies. 


II.  SOME  PRELIMINARY  CONSIDERATIONS  FOR  A  THREE- 
DIMENSIONAL  INVISCID  BLUNT  BODY  INVERSE  METHOD 


The  equations  governing  the  steady  three-dimensional  rotational  flow  of  a 
compressible  inviscid  gas  may  be  written  in  general  vector  form  as 


V  •  pV  =0 

(1) 

— >]<  — >;< 

p  (V  •  V)V  +  VP  =o 

(2) 

V  •  V(p  /p  Y)=  0 

(3) 

where  the  (  )  denotes  dimensional  quantities. 

For  a  three-dimensional  blunt  body  inverse  technique  formulated  in  terms 
of  a  curvilinear  coordinate  system,  a  coordinate  surface  corresponding  to 
the  shock  surface  is  chosen  as  a  portion  of  an  ellipse  of  rotation  about  the 
z-axis.  This  portion  is  to  lie  in  the  half- space  x  <  0,  where  the  x  =  0 
plane  extends  past  the  sonic  point  on  the  shock  surface.  The  free  stream 
velocity  vector  is  parallel  to  the  x-axis.  In  oblate  spheroidal  coordi¬ 
nates  the  metrics  are 


cosh2  | 


(cosh2  |  -  cos2  r|) 


,2  _ _ 

h2  _2 


2 

COS  T) 


C^cosh2  |  -  cos2  t|) 


,2  _ 
h3  „2 


2  r 
COS  L, 


C“  cosh2  %  COS2  T] 


(4) 
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for 


4  l  =  C  sinh  4  4  2  =  sin  t|  4  3  =  sin  4 


(4a) 


The  coordinates  4^1  4  £3  are  t'ie  curvilinear  coordinates;  values  of  4. 

r) ,  4  are  essentially  parameters  constant  on  families  of  the  4^>  £3 

coordinate  surfaces,  respectively.  Nondimensionalization  by  free-stream 
values, 


u 


p  v2 

r  00  CO 


S  =^-  (Y  -  1) 


(4b) 


is  next  performed  on  the  system  (Eqs.  1  through  3).  Note  that  the  metrics 
(Eq.  4)  are 


hr  = 


§r  +  c 


1  - 


hr  = 


1  -  t‘ 


+  2  t  +  chl  3  (ef +  o2)(i-4i) 


(4c) 


and  the  surface  4  =  const,  corresponds  to  an  ellipse  of  rotation  about  the 

z-axis;  the  surface  r|  =  const,  corresponds  to  a  one- sheet  hyperbola  of 

revolution  about  the  z-axis;  and  4  =  const,  corresponds  to  a  plane  through 

the  z-axis  (Fig.  1).  The  constant  C  corresponds  to  the  origin- to-focus 

3 

distance  of  the  coordinate  system,  or  the  radius  of  the  4  =  0  disc. 


C  is  the  common  focus  point  of  the  4  =  const.  ,  q  =  const,  surfaces;  the 
plane  z  =  0  minus  the  4  =  0  disc  is  the  surface  r|  =  0. 
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For  the  range  0  <  £  <  £  (where  £  is  taken  to  represent  a  portion  of  the 

prescribed  shock  in  the  half-  space  x<t0)0<ti<tt/2,  0  <  £  <  tt  /2,  points 

(£>r],  £)  of  the  prolate- spheroidal  coordinate  system  cover  the  space  of  the 

three-dimensional  shock  layer  (R^)  lying  between  £g  and  B(£,Tp  £),  the 

coordinates  of  a  computed  body  surface.  The  shock  £g  and  are  chosen 

such  that  the  focus  of  £g  lies  inside  of  the  computed  body  corresponding  to 

£  .  In  this  coordinate  system  (Fig.  1),  then,  Eqs.  (1),  (2),  and  (3)  become, 
s 

in  scalar  form 
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or^+Wv<Heiid=0  •  (pvi)-=o 

Vj  9(pp~v)  =  0 

8x^ 

where  v1  are  the  contravariant  velocity  components. 


(5a)- (9a) 


A  point  e(R-^nA)  is  formed  by  the  intersection  of  the  orthogonal 

(£r  t2>t3)  surfaces;  at  P,  on  the  plane  £,  =  const.  ,  unit  vectors  e^,  ^  are 
tangent  locally  to  the  £  =  const,  and  p  =  const,  traces  on  £,  =  const.  , 
respectively.  The  unit  vector  e^  lies  in  the  direction  of  X  (-L  to  the 
plane  £,). 


The  intersection  of  £,  p  lies  on  a  circle  (at  z  =  const.  )  with  the  unit  tangent 
vector  e^.  The  e^  vectors  change  direction  according  to  P.  The 

vectors  ,  V|  2>  are  1  to  4  3.  respectively;  therefore, 

e^  =  h^V^  •  W  now  dP^  is  tangent  to  the  4  ^  circle,  |  dP^  j  =  dS^  ,  then 
dP^  •  e^  =  dS^  =  dP^  •  h^V^  ,  or  dS^  =  h^d^  .  Similarly,  =  h^V^  , 

®1  =  ^1^1  ’  or  dS^  =  ^2^2  ’  =  la i i  .  The  unit  vectors  are  such 

that  =  e^  X  X  ^  ,  and  X  • 

In  this  coordinate  system  x  =  -cosh  £  cos  p  sin  t,,  y  =  C  cosh  4  cos  p  cos 
z  =  C  sinh  4  sin  p;  in  a  plane  £,  =  const.  ,  the  coordinate  surfaces  |,  p  are 
of  the  form 

2  2  2  2 

2 _  4-  _ 6 _  =  1  P _ - _  =  1  (10) 

2  2  T  2  2  1  ’  2  2  2  2  1 
C  sinh  4  C  cosh  4  C  cos  p  C  sin  p 
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2  2  2 

where  p  =  x  +  y“  ,  £,  =  0  corresponds  to  x  =  0,  and  £  =  tt/2  to  y  =  0.  In 
order  to  transform  the  (£,r|,  £,)  system  (Eqs.  5  through  9)  into  a  coordinate 
system  (<r,  t,  it)  in  which  the  shock,  body,  stagnation  streamline,  and  sonic 
surface  are  known  in  terms  of  values  of  the  new  coordinates  in  the  manner 
of  Ref.  1  (the  method  of  Ref.  1  pertains  to  axisymmetric  elliptic  regions), 
stream  functions  (Refs.  2  and  3)  ijj,  <p  are  introduced.  Note  that  the  mass 
flow  m  ,  through  a  portion  of  a  surface,  A,  covered  by  curvilinear 
coordinates,  a  and  j3,  (with  metrics,  h  and  k)  is 


=  JJ  (Vij;  X  V<f>)  ■  iTdS 

°°  An 

=  //  t(h+a®l  +  k^3®2)  *  (h*a*l  +k^0®2)]  '  “hiT  dadp 

An 

=  (4^2  '  4Ji)(9i2  '  ^ 

therefore,  for  a  plane  -x  =  const.  (AV  )  divided  into  rectangular  "boxes" 
by  z  =  const,  and  y  =  const,  planes,  -  zn_i)(yk  -  Yk_  j )  = 

(p  -  vp  “  0k  represents  the  m  in  terms  of  the  area  of  a 

particular  "box.  "  For  plane  flow  (y^  -  y^_j)  =  1,  (0k  -  i )  =  1 »  and 
(Eq.  11)  becomes  m^  =  (p^  ~  i]^)  • 

Consider  next  the  independent  variable  transformation 


I  t  \  p  /„  <r  \  _  0 _ 

\Q  sin  p/  sinh  £  sin  -q  ’  v  cos  f, /  cosh  §  cos  r|  cos  t, 

w  =  aKii^)(^r) sinhe  (12) 
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or 


Q'  JL  ,  !k 

r 

b 

L_« 

0 

°  Sfc  ' 

L  (‘-^)1/2J 

[(cz  +  lj)(i  -  6|)]1/2 

(1  ?a) 

where  a,  (3,  a',  P1  are  constants.  Through  Eq.  (12),  the  region  RjHA 

(S<  __ 

transforms  to  R^OA;  that  this  is  true  will  be  demonstrated  after  the 
following  preliminary  remarks  leading  to  a  method  based  on  Eq.  (12)  and 
generalizations  of  the  mass  flow  integral  Eq.  (11). 

For  certain  of  the  more  elementary  methods  of  the  first  type,  it  is  neces¬ 
sary  that  normal  derivatives  h^9u/9^,  hj^8v/9^,  h^9w/9^,  h^9p/9^, 
and  h119p/9^  (and  higher  derivatives)  to  the  surface  be  determined 
in  terms  of  known  tangential  derivatives  and  evaluated  from  shock  data. 
In  this  manner,  provided  a  high  accuracy  difference  scheme  is  coupled 
into  the  method,  data  is  found  on  an  adjacent  £,  surface  to  the  shock  after 
multiplying  by  a  suitable  value.  Repeating  the  evaluation  of  the 
derivatives,  the  entire  flow  field  is  then  determined,  provided  certain 

finite  difference  conditions  on  optimum  step  size,  stability,  systematic 

4 

errors,  truncation  error,  etc.  ,  are  satisfied.  These  conditions  have 


Straightforward  methods  of  this  type  for  three-dimensional  flow  will  in¬ 
variably  experience  numerical  difficulties  from  growth  of  round-off  error , 
etc.,  and  require  progressively  more  non-central  high-order  difference 
schemes  to  be  used  in  conjunction  with  the  |  set  of  derivatives  near  the 
sonic  surface,  therefore,  implying  a  much  increased  manipulation  of  num¬ 
bers  (and  techniques)  with  no  guarantee  of  convergence  to  sonic  line  values. 
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only  been  theoretically  determined  for  the  improperly  posed  Laplace 
equation  (Refs.  13  and  14);  however,  it  is  assumed  that  the  underlying 
principle  of  these  more  simple  cases  carries  over  to  the  elliptic  non¬ 
linear  system  of  equations.  A  more  satisfactory  approach  is  based  on 
certain  transformations  (Refs.  1  and  2)  initially  made  on  the  PDE  system 
in  order  to  reduce  numerical  difficulties  associated  with  the  above  method. 
In  this  method  known  sonic  line  boundary  conditions  are  naturally  intro¬ 
duced,  and,  in  addition,  flow  boundaries  and  streamlines  are  known  in  the 
transformed  system  permitting  a  reduction  in  the  order  of  the  transformed 
system.  Application  of  known  sonic  line  boundary  conditions  (M  »  1) 
are  generally  absent  from  methods  of  the  first  type.  As  both  methods 
essentially  involve  marching  methods,  clearly  generalized  considerations 
of  Refs.  13  and  14  are  of  significance.  Other  types  of  solutions  for  three- 
dimensional  subsonic  flows  are  suggested  by  Ref.  3,  in  which  systems  of 
nonlinear  ordinary  equations  for  terms  of  a  series  expansion  for  density 
and  stream  function  pairs  about  the  flow  axis  are  derived;  further  con¬ 
siderations  to  solutions  of  this  nature  will  not  be  given  in  this  report, 
however. 
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where 
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70^  afr  {h2h3s 


12  3 


afj  (hih3)  (hih2>  0 


1  /  Sh2  8hl\ 

V^ivaI7-u8^j 

1  /  ah3  8hl\ 


1  /  9hl  9h2\  1  /  9hl  8h3\ 

e?t2^w2-v^)  0 

0  tttt-  (v  -T—  -  w  -t—  ]  0 

h2h3  \  8^3  9^ 2 / 


1  (  8h3  8h2\ 

0 


and  boundary  conditions  on  the  prescribed  =  const,  shock  are,  for  a 
perfect  gas  at  constant  Y,  at  a  point  P(|  ,  p,  C),  for  V  =  7!  e  h  f11 


i  u 

‘‘■sf  ' 


cos  r|  sin  f,  sinh  £  g  (Y  -  1 ^  +  2 

2  *  i  ■  „ 2  __  \1  / 2  ,,,  ,  i  »»i2 


1  h^|sinh  4g  +  sin  r| )  (Y  +  1)M^ 


(10a) 


P=  '  [>  +  7TT  K>,  N  -  ')] 


O  w  V 

S  COS  L  /ii  it\ 

f  =  r—  = - tt  “  11a- 12a) 

h3  h3 


9  v  sin  t)  sin  £  cosh  £  . 

_  S  SI 

1  -  h"7  /  .  “  2  :  7  •  2  iJ72  h7 

2  (sinh  §  +  sin  r| ) 


p 

Ps'lv-  DM^Nt2 


(1 3a- 14a) 
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and  the 


where  ^  =  M^jcos2  r|  sin^  £  sinh^  |g/(sinh^  |g  +  sin^  q) 

corresponding  entropy  is 


S 

s 


-  S 


CD 


log 


1 

r  2  H 

-V 

[i  +  2Y  /  m2  \1 

<^  +  1Ko,n 

Y+l 

L  \  oo.  n  )J 

.'V-  1)<N  +  2 

(13) 


while  on  the  body,  b,  S  =  S  I  =  const. 

1  'max. 


(13a) 


d£^cosh^  |  -  cos^  q 

1/2- 

=  C 

b 

\*Z 

t  2 
cosh  5  -  cos  T) 

1/2-1 

u  cosh  £ 

_ 

_ 

V  COS  T) 

J 

d£, ^  cosh  £  cos  t] 


w  cos  £ 


(14) 


Suppose  now  that  Eqs.  (5)  through  (9)  are  subject  to  independent  variable 
transformations  (Eq.  12)  by  means  of  two-component  stream  functions, 
but  not  put  into  Eqs.  (5a)  through  (9a)  form,  as  the  calculation  scheme 
proposed  does  not  require  this  form.  Then  conditions  on  the  sonic  surface 
may  be  determined  by  characteristics  methods  so  that  RPlA  conditions  are 
known  on  a  closed  boundary. 

In  that  curvilinear  coordinates  are  used  throughout,  differentiations  with 

respect  to  coordinates  are  best  done  expressing  vectors  in  their  contra- 

3 

variant  form  (e.  g.  ,  V  =  ^h^f n,  where  f n  are  components  of  a 
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contravariant  vector);  this  provides  a  systematic  formalism  for  the 
effects  of  differentiation  of  the  "unit  vectors.  "  In  general,  then,  for 
the  velocity, 


9V  _  ^  A  , 

3F  -  A  6nhn 
n=l 


9f 


3 

+  E  fr 

m=l 


l|  nl 

lm  jl 


(14b) 


or  the  components  of  the  contravariant  vector  which  correspond  to  the 


derivative  of  the  ordinary  vector  with  respect  to  are 
3  J 


f  T  ;  =  +  E  symbo1  {m^j  1  benotes  tbe  Christoffel 


symbol  of  the  second  kind,  generally  formed  from  those  of  the  first  kind 

c  |  k  \  kar..  u  f  i  i  1  /^ik  ^®jk  ^ij\  j  ka 

from  |  =  g  [ij,  a]  ,  where  [ij,  k]  =  -^1 - r-  +  — *4 - ^1 ,  and  g 

J  \  9xJ  9x  9x  / 

is  the  contravariant  tensor  g1^  =  G1!/  |g^  |,  GXJ  is  the  cofactor  of  the 


element  g. .  in  g. 

°1  1  1  °l 


ij 


3iJ 


Provided  transformations  are  next  taken  with  the 


covariant  law,  then  (for  scalars),  with  Eq.  (12a), 


3(  ) 

9(j-  9t  9tt 

9(  ) 

^7  ^7  ^ 

9  <r 

3(  ) 

9  or  9t  9tt 

9(  ) 

^ 

9(  ) 

9cr  9  t  9tt 

9  (  ) 

**7 

jt| 

CD 
l TP 

1 

1  9tt 

(15) 


where  the  column  vectors  of  J  represent  covariant  vectors.  In  terms  of 
these  quantities,  then,  the  derivatives  of  the  boundary  conditions 
(Eqs.  10a  through  14a,  and  13)  may  be  expressed  in  terms  of  the 
Christoffel  symbols,  using  Eqs.  (4a),  (4c),  and  (14b). 
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At  the  sonic  surface 


holds.  Also,  from  Eq.  (12a)  at  the  shock 


(17) 


The  constants  K^,  are  defined  at  the  shock  by  Eqs.  (1 8)  and  ( 1 9);  at  the 
shock  front  in  either  coordinate  system,  (cr,  t,  tt)  or  (l^,^,^),  f  is  known 
from  Eq.  ( 1  2a)  and  the  related  inverse  relations.  The  relation  Eq.  (17) 
holds  throughout  (R^HA)  and,  when  used  numerically  with  Eq.  (16),  yields 
a  method  whereby  the  sonic  surface  may  be  computed  (Refs.  1  and  3). 
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To  describe  the  transformation  (Eq.  12),  consider  the  quantity  |e  |  as  a 
small  positive  constant,  |«  j  «  1.  The  quantities  a,  (3,  a',  (31  are  also 
small  constants,  or  scale  factors,  on  t,  it  (Eq..  12  and  12a).  By  con¬ 
struction  of  a  suitable  approximation  surface  (to  be  described  later  in  this 
section)  through  the  neighborhood  of  z  =  0,  y  =  0  on  the  -x  =  const,  plane 
(stagnation  point),  and  on  the  stagnation  streamline  though  this  point  to 

the  =  const,  shock,  ijj  and  0  are  of  the  order  of  |  e  |  on  this  surface. 

1  s 

On  the  body  surface,  where  iji  =  0  and  0  =  0  simultaneously,  the  i|j,  <p  are 
of  the  order  of  |  e  |  on  the  approximation  surface  Q  in  the  neighborhood  of 
any  point  on  the  body  surface.  Away  from  these  regions,  >>  |«  |, 

0  >>  |«  |  and  may  be  considered  simply  as  i)j,  0  respectively.  On  the 
shock,  from  Eq.(12),  assume  (ar/sinr|)  =  const.,  (pcr/cos  !!,)  =  const, 
so  that  in  terms  of  (<r,  r,  it)  coordinates,  it  =  (const.  )  sinh  |.  Then,  as 

=  const.  ,  the  prescribed  three-dimensional  shock  in  oblate- 
s 

spheroidal  coordinates  transforms  to  it  =  (const.  )  in  the  new  system 
(Fig.  2).  As  0  <  q  <  it/2,  0  <'  t  <  const.  ;  0  <  £  <  tt/2,  Ocacconst., 
outside  17.  On  £  ,  coordinates  (C, ,  -q )  determine  a  point;  corresponding 
to  this  point  is  (t,  <t)  on  it  =  (const.  ).  For  a  locally  nonvanishing  Jacobian, 
then,  a  portion  of  the  shock  A(£g,q,  £)  may  be  mapped  into  a  region 
A'(tt,  (t,  t)  in  the  new  coordinate  system.  Provided  next  the  box-like 
regions  (see  Eq.  11  for  simplified  discussion)  on  the  -x  =  const,  plane 
are  projected  onto  the  £g  =  const,  shock  (Fig.  1),  then  on  the  shock  surface 


~  sinh  ^  sin  q  ~  z 


0~  cosh  |  cos  q  cos  t,  ~  y 


(18) 

(19) 
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r^n  fz  n  f9  k  fy  k 

or  /.  dj/  ~J  dz  and  d0  —  /  dy  for  regions  away  from  the 

•^n-l  ^  Zn- 1  ^k-1  yk-l 

f2  trace  on  |g.  The  projected  box-like  regions  do  not,  in  general,  conform 

to  the  curvilinear  r|,  £,  system  on  the  shock  surface,  4S-  Immediately 
downstream  of  the  4S  shock,  relations  Eqs.(18)  and  (19)  are  no  longer 
valid.  For  a  streamline  of  the  flow  field  (which  lies  along  the  inter¬ 
section  of  Jj  =  const.  ,  <p  =  const,  surfaces  on  a  common  value  entropy 
(S/R)  line  imbedded  in  both  surfaces),  the  conditions  become  i|/  =  const.  , 

0  =  const,  on  a  general  streamline,  or,  from  Eq.  (12), 


■h 


CLT  = 


const. 


const. 


s 


mhl; 


^  cosh  4  cos  T| 


(20) 


and  in  the  neighborhood  of  the  body 


_  e  /2 

aT  sinh  4 


l«  1/2 

9<T  cosh  4  cos  q 


(21) 


as  O(iJj)  —  0(  |  e  [ );  0(0)  ~  0(  |«  | ),  near  the  part  of  12  over  the  body. 


In  that  the  stagnation  streamline  and  body  stream  surface  have  an  entropy 
(S/R)g,  it  is  necessary  that  a  suitable  representation  be  found  for  com¬ 
binations  of  the  new  coordinates  to  represent  the  adjoining  stagnation 

streamline  and  body  as  a  surface  with  similar  properties;  to  facilitate  this 

I  l 2 

representation,  a  surface  consisting  of  a  tubular  region  of  radius  |e  | 
surrounding  the  stagnation  streamline  joining  a  surface  that  envelopes 


^Note  that  the  regions  of  12  are  chosen  such  that  the  |e  |  inequalities  are 
satisfied  for  a  definite,  given  |«  |  value. 
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the  body  (within  j  e  j  of  the  body)  is  assumed  to  represent  the  (S/R) 

5 

surface.  In  this  manner  a  continuous  approximate  surface  is  con¬ 
structed  that  lies  arbitrarily  close  to  the  actual  stagnation  streamline 
and  body  surface,  by  suitable  choice  of  an  |*  |  value.  In  the  limit  as 
j  (  ]  -*  0,  the  actual  stagnation  streamline  and  body  surface  are  exactly 
approached.  Now  in  the  region  of  the  stagnation  streamline, 

O  j  cos  £,  |  ~  O(q)  —  0{  j  e  j ),  0  <£<£g  (as  the  focus  of  the  coordinate 
system  is  assumed  to  be  within  the  body)  so  that  in  the  limit  lei-*-  0, 
r  -*  0 ,  <r  -*  0,  and  u  appears  of  the  order  of  magnitude  aprnax  I  sinh  |  | 

U  -  4S  i  <r  6 

in  the  region  of  the  stagnation  streamline  just  behind  the  shock  (6  is  a 
small  positive  constant).  This,  then,  in  (<r,  t,  it)  coordinates,  appears 
to  be  the  it  axis,  tr  <  (const.  )max  |  sinh  §  |  .  In  the  vicinity  of  the  body 

lf-esl<6- 

surface,  §  >  0,  ir/2  >  q  >  0,  tt/2  <  £  <  0,  away  from  the  tubular  region  and 
near  £2.  The  transformation  system  near  £2  is  then, 


/  T  \  sinh  £  sin  q  = 
\  sin  r\)  5  1 


(  \ 
\cos  f,  ) 


cosh  |  cos  q  cos  £  =  |  e 


(22) 


(23) 


IT 


(^2L\(_§E_)  sinh  4 

\sin  q/\cos  f,  / 


(24) 


2 

Suppose  that  |  is  negligible  with  respect  to  |  €  I  and  1  e  |  is  taken  to  be 
a  small  positive  constant  in  the  neighborhood  of  the  body.  From 
Eqs,  (22),  (23),  and(24),  then,  0(ar/sin  q)  ~  0(  |  e  |  /m),  0(  a/cos  Q  ~ 

0(  |  <  |  /n),  where  n,  m  may  be  chosen  >1  for  large  £>  (3  and  for  conditions 
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sufficiently  far  from  the  stagnation  point  region.  Then 

i  1 2 

O(ir)  ~  0[(  |  «  |-/  mn)  sinh  £]  .  On  the  approximate  surface,  then,  for 
conditions  applicable  near  the  forward  part  of  the  stagnation  streamline, 
part  of  the  it  axis  is  reproduced;  further,  for  conditions  on  the  body 
surface  sufficiently  far  from  the  intersection  point  (p  >  0,  t,  >  0,  £  large) 
the  it  surface  becomes  asymptotically  close  to  the  tt  =  0  plane.  The 
region  of  0  in  the  vicinity  of  the  stagnation  point,  or  the  juncture  of  the  two 
approximate  surfaces  is  yet  to  be  examined.  Initially,  consider  the  z,  y 
planes  tangent  to  the  outside  of  the  |e  |  tube  and  those  that  describe  a 
square  within  the  tube  on  the  plane  -x  =  const.  ;  the  ijj ,  0  planes  lying 
between  these  surfaces  (in  this  approximation)  then  generate  an  infinity 
of  (ij;,  0)  pairs  such  that  any  vjj ,  0  pair  in  this  region  forms  streamlines 
lying  both  on  the  |«  |  tube  and  the  approximation  stream  surface  within 
|  €  |  of  the  body  surface  while  still  imbedded  in  ft.  On  the  approximation 
body  surface,  then,  for  a  \e  j  ^  0,  specifying  a  point  (£,r|,  £)  leads  to  a 
point  (cr,  t,  tt  «  0)  in  the  new  coordinates  for  regions  sufficiently  far  from 
the  stagnation  point.  In  the  region  near  the  stagnation  point  on  ft,  suppose 
that  0(£)  ~  0(  I  e  | ),  then^  O(tt)  ~  0(a(3  |  e  | ).  Here  it  is  assumed  that  in  this 
region  the  0(sin  0  )  ~  0(  |  e  | ),  0(cos  ir/2  )  ~  0(  |  e  | ).  The  value  of 
sinh  |  for  £  »  1  is  effectively  of  the  order  of  1/2  e^;  for  £  small,  then, 

0( s inh  £)  ~0(-l/2  e  ^),  (£  >  0).  For  (ar/sin  r|)  and  (p^/cos  £)  of  the 


Here  it  is  assumed  that  the  coordinate  point  (0,  ££>  £^)  lies  in  the  region 
of  the  stagnation  point;  note  also  that  the  stagnation  point  lies  interior 
to  ft  in  the  neighborhood  of  that  point. 
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order  of  J  t  \  t  0  for  r|  >  0,  £  <  tt/2  on  the  body  (Fig.  1),  then 

O(tt)  —  0(  [ «  |  ).  The  constants  a.  p  may  be  suitably  chosen  according  to 

the  degree  of  magnification  of  the  cr ,  t,  tt  coordinates  desired  for  a 

7 

particular  problem.  For  tjj,  <p  increasing,  then,  the  value  of  it  is  >0, 
consistent  with  the  streamlines  passing  between  the  shock  and  body  in  the 
(<r,  t,  tt)  space.  These  streamlines  originate  at  intersections  of  ijj  =  const.  , 
<p  =  const,  surfaces  on  an  x  <  0  plane,  or  equivalently,  the  mappings  of 
these  intersection  points  onto  the  tt  =  const,  plane  in  the  (cr,  t,  it)  space. 

For  £  large,  then,  sinh  |  is  effectively  larger;  therefore,  the  scale  of 
the  tt  coordinate  is  exaggerated.  The  surface  corresponding  to  £  =  0  is 
the  surface  of  disc  r  =  C;  by  requiring  the  stagnation  point  to  be  close  to 
the  focus  of  the  coordinate  system,  then,  the  three  parts  of  the  approxi¬ 
mation  surfaces  satisfy  the  properties  previously  designated.  In  this 
manner,  the  forward  part  of  the  stagnation  streamline  just  past  the  shock 
transforms  into  part  of  the  Tr-axis;  the  part  of  the  stagnation  streamline 
closer  to  the  stagnation  point  lies  on  a  surface  O(tt)  ~  0(  [c  |),  and  the 
surface  covering  the  body  within  je  |  possesses  the  property  that 
O(tt)  ~  0(  j  6  |).  In  the  quadrant  of  the  (o',  t,  tt)  space  into  which  the 
quadrant  of  the  physical  (£,q,  £,)  shock  layer  space  has  been  mapped,  this 
surface  is  then  of  the  general  nature  of  a  hyperbola  of  revolution  with  the 
it,  n-,  and  t  axes  as  asymptotes.  Note  that  from  Eqs.  (16)and(17),  and  the 


Obviously  the  magnitudes  of  these  constants  are  determined  from  geo¬ 
metrical  considerations  based  on  the  chosen  £  shock,  M  ,  C,  etc. 

s  oo 
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results  of  Section  IV,  the  sonic  surface  may  be  constructed;  the  closed 

•  ■V 

region  R^HA  therefore  consists  of  the  given  shock,  known  sonic  surface, 
and  Cl  in  (g-,  t,  it)  space  (Fig.  2).  From  Section  IV  it  is  also  possible  to 
construct  stream  surfaces,  a  priori,  through  this  region.  Note,  also, 
that  for  transformation  of  vector  quantities,  Eq..  (14b)  must  necessarily 
modify  the  form  of  Eq.  (15)  due  to  the  presence  of  contributions  from  the 
differentiated  unit  vectors.  This,  then,  roughly  describes  the  bounding 
surfaces  of  the  physical  shock  layer  in  terms  of  surfaces  of  the  (o-,  t,  it) 
space,  provided  certain  assumptions  are  made  concerning  approximation 
surfaces  to  the  actual  (S/R)^  surfaces.  The  Jacobian  of  the  transformation 
(for  scalars)  is 


3c r  3cr  3c r 

aiy  sf; 

9t  3t  3t 

5 XI  w> 

I  3tt  3  it  3tr 

and  J  and  the  inverse  j  =  J  ^  are  to  be  values  bounded  away  from  0,  co  in 
the  region  bounded  by  the  shock,  sonic  surface,  and  the  continuous  surface 
enveloping  the  stagnation  streamline  and  body,  in  order  to  guarantee  a 
continuous  one-to-one  mapping  of  points  in  one  space  to  the  other.  This 
condition  may  require  adjustment  of  the  arbitrary  constant  values  intro¬ 
duced  throughout  this  Section. 


=  J  ;  cr,  t,  ne  C 1 


(25) 
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III.  THREE-DIMENSIONAL  STREAM  FUNCTIONS 
OF  CLEBSCH  AND  MORE  GENERAL  TYPES 


For  the  inviscid  three-dimensional  rotational  flow  of  a  compressible  gas 
in  a  oblate  spheroidal  coordinate  system,  the  streamline  equations  are 
obtained  as  integral  curves  of  the  energy  equation 


u 


cosh  £ 

3 

/  _E_\ 

(  2  \1/2 
Icosh  £  -  COS  Tj  f 

■ 

lp\) 

cos 


T7T 


^cosh^£  -  cos^qj 


+ 


w 

C 


COS  L, 

cosh  £  cos  T)  3 


(5) 


0  (26) 


3  'p 


as 


,1/2 


d£^cosh^£  -  cos^qj  d£.,^cosh2£  -  cos^qj 

- - — - -  -  C  ~ 


u  cosh  £ 


Zt  2 

i  t,  -  co 

V  COS  T| 


1/2 


=  C 


d§ ^  cosh  £  cos  q 
w  cos  £, 


(27) 


The  general  integral  of  Eq.  (26)  is  of  the  form  p/ p^  =  f[4i(£,  q,  u,  v,  w), 

8  v 

0(£,  q,  £,,  u,  v,  w)]  where  f  is  determined  from  initial  values  of  p/pY  on  the 
streamline  intersection  of  ij;  =  const.  ,  <p  =  const,  surfaces  with  a  shock, 


See  Eq. 


(17). 
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here  taken  as  a  portion  of  the  coordinate  surface  £  =  const,  for  x  <  0.  In 
that  both  tjj  and  <f>  satisfy  Eq.  (26)  individually,  it  follows  that 


cosh  £  i  , 

u - 2 - +  v 


cos 


7  z  r?71  ^i  I  z  r?71  ^2 

|cosh  £  -  cos  r)j  C^cosh  £  -  cos  r|j 


w  cos  f,  .  _  - 

C  cosh  £  cos  r| 


(28) 


and 


u 


cosh  g 


2  i  TT71^ 

^cosh  £  -  cos  r)j  C^cosh  £;  -  cos  qj 


cos  T| 


w  cos  £  _n 

C  cosh  £  cos  t) 


(29) 


Solutions  u,  v,w  of  Eqs.  (28)  and  (29)  are  determined  to  within  a  propor¬ 
tionality  factor:  by  requiring  these  values  to  satisfy  the  equation  of 
continuity 


af^P^cost^)  (cosh2£  '  cosS)  ] 

+  a!l[pv(§^)  (c°sh2!  -  cos2pj  7j 
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9 

it  follows  that  u,v,w  are  determined  in  terms  of  the  tjb  <P  stream 
function  pair  as 


1  cos  r|  cosh  £ 

P  C(cosh^£  -  cos^r)) 


-1 


(%%  -  \\) 


(31c) 


which  are  the  components  of  V  =  1  / ptV+xV^).  In  order  to  transform  the 
basic  quasi-linear  first-order  system  (Section  I)  into  another  first-order 
system  by  means  of  dependent  variable  transformations  based  on  Eq.  (31), 
it  is  necessary  that  inverse  relations  of  the  form 


^1  =  l^1^l,^2’^3’U’V’W'P* 


(£r  i2>  u,  v,  w,  p)  , 

2  2 


3 


3^r  ^ 


3’ 


u,  v, w,  p) 


=  ^1(^1’^2,^3,U,V,W,p)  ’ 


^  =  "^g  (cr ^3. u- v>w» p)  *  0|  =  i>g  (€1.e2.e3.u,v,w>P) 


Q  2 _ 

7  Alternatively,  V<p  -  V^x(  *  p/ |v4*|  )V;  note  that  if  4^  is  given,  then 


-I 


(VijJxpV)  •  dx 

.  Vvp 
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are  known.  Clearly  the  system  (Eq.  31)  is  not  sufficient  for  this  purpose 

as  implicit  function  theory  yields  only  relations  (for  example)  of  the  form 

cos  r|94u  /9(pu)  -  C4J*  cosh  /9(pv)  =  0,  which  lead  to  general 

*2  ’3  M  53 

integrals  only  reasserting  conditions  Eqs.  (28)  and  (29).  Additional 

assumptions  are,  therefore,  necessary  for  determination  of  suitable 

inverse  relations  from  Eqs.  (31a),  (31b),  and  (31c),  provided  this  system 

is  to  be  used  to  reduce  Eqs.  (5)  through  (9). 

Transformations  of  the  type  V  =  l/p(V^xV0)  are  Clebsch  transformations. 
An  alternative  type  is 


puh2h3  =  Ml  (32) 

pvhih3  =  iMi +  U^)  (33) 

Pwhlh2=f|^  (34) 

where  the  metrics  are  as  given  in  Eq.  (4c).  Provided  pV  =  curl  A  then 
the  existence  of  real  ip,  <p  surfaces  in(ROA)  becomes  the  problem  of  deter¬ 
mining  real  values^  for 


^For  aVxB  =  VipxV^*,  B.2h.  =  >  B?2h?  =  ip0t  -  \it  , 

11  *1  51  ^ 

B  2h^.  =  -  0ipt  >  illustrating  the  difficulty  associated  with  a  form 

■5  -3  53  S3 

corresponding  to  Eqs.  (35)  and  (36)  for  this  transformation.  The  forms 
Eqs.  (35)  and  (36)  are  nonunique,  however,  as  A.  terms  are  determinable 
to  within  an  arbitrary  gradient. 
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^  h3A3  ~f  Si 3  (h2A2)  d^2 

e1.e2.e3*(R1nA) 

0=  -hiAi+/af^(h2A2)d^ 


(35) 

(36) 


For  the  ip.  <p  stream  functions  of  Eqs.  (32)  through  (34)  to  satisfy  Eq.  (26) 
individually,  as  (Eqs.  28  and  29),  then 


"Vlj  h2%  h3^3  uh3 

hl%  =  H2%  =  h3^3  = 


(37a, b,c) 


from  Eqs.  (32)  through  (34).  This  follows  from  the  requirement  that  a 
point  on  a  streamline  is  imbedded  on  both  <f>  and  iJj  surfaces. 

On  the  4*  =  const,  sheet  and  the  <j>  -  const,  sheet,  respectively,  additional 
relations  follow  from 


d^lhl  +  h21^  d^2h2  +  d^3h3  =  0 

1  2  3 

blS^^l  +  h2l0|2de2h2  +  h31^3d|3h3  =  0 

So  that  with  Eqs.  (28)  and  (29), 


(38) 

(39) 


hihll*z1 


2  2 
-h^dl^vh^h^  +  h^d^wh^h^ 

2  2  bi^3 

h^d^vh^h^  -  h^d^uh^h^ 


(40) 
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(41) 


h2h3S2 

2  2 

-h^d^wh^h^  +  h^d^uh^h^ 

h 

53 

hld^lvh2h3  '  h2d^2hlh3 

-hjd^vhjh^  +  ^12d^2w^1l^13 

?3 

2  2 
hjd£  lv^12^13  "  ^2d^2u^^3 

h2h3^e2 

2  2 

-h^d^-wl^h^  +  h^d^h^h^u 

h 

53 

hld^lvh2h3  "  h2d^2hlh3U 

(42) 


(43) 


on  these  constant  valued  surfaces.  Also,  for  orthogonality  of  i|;,  <t> 
surfaces, 


.  V«2 .  %%  „ 

T 2—  +-TI—  +— =  0 


n 


1 


(44) 


which  generally  is  satisfied  only  in  the  free  stream  and  on  the  shock,  but 
is  probably  too  restrictive  in  and,  therefore,  will  not  be  assumed 
(Fig.  1).  Then  it  follows  that  in  R^,  outside  of  fi, 


-28- 


% 

% 

% 

% 

<pt 

<Pp 
=  3 

_  _j 

(uh3)(pwh1h2)^  +  (puh2h3)^  (wh^ 


(uh3)(pwh1h2)^  -  (puh2h3)^  (wh^  puh2h3^-1^ 

puh2h3 


(uh3)(pwh1h2)^  -  (puh2h3)^  (wh^ 


(uh3)(pwh1h2)^  -  (wh1  )(puh2h3)4  Puh2h3 


(uh3)(pwh1h2)^  +  (puh2h3)^  (wb^) 

( u.h3 )( p wh  1  h2 ) ^  -  (puh2h3)^  (wh^)  pwhlh2 

pwh.^h2 

(uh3)(pwh1h2)^  +  (puh2h3)^  (wh^ 

_pvhlh3  +  puh2h3  (uh3)(pwh1h2)t  -  (puh^k  (wh^ 


(45) 


At  the  shock,  these  relations  are  explicitly  given  for  the  4b  <p  derivative 
column  vectors. 


Note  that  in  Ref.  1,  the  expression  Eq.  (37)  for  the  axisymmetric  case  is 
unnecessary,  as  -  i^^k/R  =  0;  for  the  three-dimensional  case, 

however,  Eq.  (37)  is  an  additional  relation  that  must  be  satisfied. 

Note  that  in  Eq.  (45)  the  evaluation  of  derivatives  is  to  be  accomplished 
by  expressing  the  velocities  in  their  contravariant  form,  then  operating 
with  Eq.  (14a),  in  part,  for  the  complete  derivative.  This  process  has 
been  incorporated  into  Eq.  (45). 
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Note  on  the  inverse  relations: 


The  expressions  in  Eq.  (45)  contain  certain  derivative  terms.  Now 
consider  n  sets  of  values  given  a  priori  u  v  ,  w  ,  §  ^  ,  §1,  <  63  > 

etc.,  in  Eq.  (45)  with  metrics  known  and 
derivatives  approximated  by  some  discrete  approximate  expression 
involving  (1),  .  .  .  ,  (n)  quantities  and  the  appropriate  metrics.  Then  values 
v^2\  .  .  .  ,  ir^n^  may  be  calculated  for  the  first  set  of  values,  respec¬ 
tively,  and  also  .  .  .  ,  (r^n^,  ,  In  terms  of  these  quantities, 

the  partials  may  be  represented  accurately,  locally,  in  the  vicinity  of  a 
particular  set  of  values  by  this  "finite  difference"  type  approach.  In  a 
certain  small  region,  then,  derivatives  have  a  certain  "best  approxi¬ 
mation,  "  as  a  number  of  points  have  been  used  for  their  determination. 


Now  for  two  particular  sets  of  values,  (j),  (j+1),  with  the  best  derivative 

approximation,  a  box  bounded  by  planes  tt^  =  const.  ,  =  const.  , 

<r^  =  const.  ,  <r^  +  ^  =  const.  ,  t^,  +  in  R  exists.  In  the  open  region 

(k)  (k)  (k) 

N  (for  which  these  planes  are  the  closure),  then,  the  u'  ,  v'  ,  w'  '  eN 
are  to  be  considered  as  parameters  for  the  inversion  of  Eqs.  (62) 
through  (65)  for  a  particular  set  of  u^,  v^,  w^  values  1  £  t  —  n.  If 
next  u(t),  w^  are  varied  and  the  process  repeated,  then  the 

dependence  upon  the  parameters  can  be  ascertained  by  the  variances  in 
surface  shape,  average  values  within  the  box,  etc.  With  this  information 
the  inverse  relations  suggested  later  can  be  performed  numerically  with 
respect  to  the  u,  v,w,  p,  parameters.  See  also  Eqs.  (62  through  65). 
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Consider  the  special  case  of  plane  flow  in  rectangular  coordinates  at 
constant  density  (w  =  0,  wz  =  0,  p  =  1 ,  hj ,  h 2>  h3  =  1,  =  °)>  then 


— 1 

-e - 

X 

J 

-V 

4* 

y 

u 

0 

<p 

0 

K 

0 

y 

i 

-o- 

1  N 

0 

(46) 


Note  also  that  for  p  =  1,  h^  =  1,  the  irrotationality  conditions  for  the 

2  2 

stream  function  pair  are  of  the  form  V  ^  0;  V  i|)  -  0  with  <f>  =  ijj  . 

X  z 


With  these  relations,  transformation  from  the  curvilinear 

space  to  <r,  t,  it  coordinates  by  means  of  Eq.  (15)  is  now  possible.  Now 

the  metrics  of  the  two  systems  are  related  by 


,  2  ,  2  ,t2  ,  ,  2  ,t2  ,2  ,t2 

ds  =  hjd|j  +  +  hjd|3 


=  h^dx2  +  h2d<r2  +  h2dir2  (48) 
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and  for  orthogonality  of  the  two  coordinate  systems, 


a^i 

9^1 

I ae2 

*2  , 

343 

ae3 

9<r 

9t 

9cr 

9t  + 

9  cr 

8t 

9g1 

9^1 

,  ^2 

H2 

9^3 

ae3 

9<r 

97T 

8<r 

9ir 

9c 

9tt 

9^1 

^1 

,  *2 

ae3 

8^3 

=  0 


=  0 


9t  9tt  9t  9tt  9t  9it 


(49) 


(50) 


(51) 


and 


’2  Ubl 
9  cr  9c r 


9t  9t  ,  9tt  9tt 

9e25f7  8^35, 

=  0 

(49a) 

9t  9t  ,  9  it  9  it 

3£2  a^3  3^  3^3 

=  0 

(50a) 

9t  9t  ,  9tt  9tt 

9gj  8g3 

=  0 

(51a) 

The  elements  of  the  J  matrix  (Eqs.  15  and  25)  then  become,  with  Eqs. 
(45)  and  (12A), 


9<r 


"■Si 


P'  1pwh1h2 


1 


(c  +  b>  [(i  -  e|)(c*  +  tJ) 


172 


(uh3)(pwh1h2)^  +  (puh2h3)^  (whj) 

(uh3)(pwh1h2)^  -  (puh2h3)^  (wh^ 

2  2 


(52) 
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(53) 


j£_  _  °^2 

a£- 


P,'1pwh1h2 


f1-^)  [(c2  +  e^)(i  -  e|) 


172 


9<r 

all 


-P 


_(c2  +  sf)(  1  -  s^)]I/z 


X 


(uh3)(pwh1h2)^  +  (puh2h3)g  (whj) 

pvh1h3  puh2h3  -  (puh2h3)^  (wh^ 


(54) 


3t  _  1 

517-61 


(uh3)(pwh1h2)^  +  (puh2h3)^  (whj) 

a  puh2h3  (uh3)(pWh1h2)^  -  (puh^)^  (whj) 


(55) 


9t 

Hr, 


puh2h3 


C56) 


(_j  (uh3)(pwh1h2)^  -  (puh2h3)^  (wh^ 

9^3  ”  T7~  puh2h3  (uh3)(pwh1h2)^  -  (whj  Xpuh^)^ 


(57) 


9  IT  _  TT  IT  9t  ,  IT  9<r 

Hl~il  r9e2  0-9^ 


(58) 


9lT  _  __TV_  TT  9t  ,  TT  9  (T 

9 ^2  ^2  T  9^2  17  9^2 


(59) 


9lT  _  _TT_  TT  9t  ,  TT  9tT 

9£3  ■  *  |3  T  9£3  <t  9^3 


(60) 
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in  which  coupling  exists  between  the  (£j,|2,£3)  and  (ct.t.'tt)  systems. 
Note  that  the  order  of  the  system  has  been  preserved  as  Eq.  (45)  is 
known,  and  the  Eqs.  (52)  through  (60)  are  i|j,  <p  independent.  Now  the 
space  R  (locally  Euclidean)  formed  by  £2,~£3  *s  to  maP  (locally)  onto 
the  totally  Euclidean  space  R  (<r,T,ir)  in  Fig.  2.  For  the  totally 
Euclidean^  space  of  (ir.T.ir),  it  is  assumed 

~2  - 1  ~2  - 1 

h^  =  =  const.  >0  ,  h^  =  =  const.  >  0  , 

—  l 

h^  =  =  const.  >  0  (61) 

so  that  all  components  of  the  curvature  tensor  are  zero  throughout  R. 
Then,  from  Eqs.  (48)  and  (52)  through  (60), 

(uh3)(pwh1h2)^  +  (puh2h3)^  (wh^ 

r  a  (puh2h3)  (uh3)(pwh1h2)^  -  (puh^Jg  (wh^  ^  ^ 


± 


a' _2(puh2h3)2 


X  1  + 


(uh3)(pwh1h2)^  -  (puh2h3)^  (wh^ 
(uh3)(pwh1h2)^  -  (wh1)(puh2h3)^ 


-.2 


1/2 


(62) 


11 


All  the  Christoffel  symbols  vanish  for  constant  metrics. 
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1 


here  is  chosen  such  that 


K^a'2  >  S^Puh^ 

i 

throughout  R,  and  (arbitrarily),  t  >  0  in  the  mapped  octant.  The 
constants  here  also  are  to  be  such  that  Eq.  (25)  conditions  are  met. 


■(uh3)(pwh1h2)^  -  (puh2h3)^  (wh^-j 

1  +'  (uh3)(pwh1h2)^  -  (wh1)(puh2h3)^ 


Similarly, 


P'  ^puh^hj 


[(i  -  «z)(c2  +  €?>] 


172 


g  (uh3)(pwh1h2)^  +(puh2hj)g  (whj)  g 

1  -  £2  (uh3)(Pwh1h2)^2  '  (puh2h3)g^(whj)  c2  +  g2 


f  6? 

*5  1 

-i 

± 

\  ,  <\  1 

_(c2  +  |; 

if  (i  -  ei)  J 

>2  +  g^  (i-52)2. 

P'-2(puh2h3)Z  [  |2  (uh3)(Pwh1h2)|i  *(puh2h3)Ci(whl) 

+ 


(c2  +  |‘)(i  -  I l) 


_2 


,1/2 


2  .  ,2 


1  -  |2  (uh3)(pwhjh2)g  -  (puh2h3)g  (whj)  C  +  gj 


-  sr 


(63) 


where 


L(c2  +  ^i) 


P'-Z(puh2h3)2 


(i-#(c2  +  e?) 


(uh3)(pwh1h2)g  +(puh2h3)g  (whj) 
(uh3)(pwh1h2)g  -  (puh2h3)g  (whj) 


P'-2(puh2h3)2 

+  (c2 ;  ^  j)(i  -  if) ' *2 


(c2  +  g2)(i  -  § 


pvhlh3  ~  puh2h3 


(uh3)(pwh1h2)g  +  (puh2h3)g  (whj) 


(uh3)(pwh1h2)g  -  (puh2h3)g  (whj) 


\2 


(64) 
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Again,  is  chosen  such  that  cr  is  real  in  R,  and  cr  >  0  in  the  transformed 
octant,  representing  the  physical  octant  in  R.  In  the  same  manner  as  the 
above,  with  (51a) 


where  the  above  partials  are  to  be  evaluated  from  Eqs.  (62)  through  (64), 

and  (52)  through  (57)  expressions.  For  values  of  u,  v,  w,  then 

t  =  t(£  1(  £2,  £3;  u,  v,  w),  cr  =  £3;  u,  v,  w)  and  it  =  tt(£  l ,  £3  ;u,  v,  w); 

2 

provided,  for  fixed  u,  v,  w  at  a  point,  J  ^  0,  j  ^  0  and  or,  t,  rreC  ,  and  cr,  t,  tt 
are  well  defined  for  all  £2>  £3«Rj,  then  an  inverse  exists  for  fixed 
u,v,w,  i.  e.  ,  ^  =  |^(o-,  t,  irju,  v,  w),  £2  =  £2(o-,t,it;u,v,w),  and 
^  =  ^3(cr,  r ,  tt;u,  v,  w).  This  procedure  was  suggested  earlier  in  Section 
III;  note  also  that  p  may  be  included  in  the  same  category  as  u,  v,  w,  in 
this  inversion. 

The  constants  K^.K^K^  change  scale  sufficiently  in  the  (ct,t,tt)  system 
such  that  real  values  of  these  coordinates  are  defined  for  real  values  of 
^l’^2’^3  coorc^nates'  The  complexity  of  the  suggested  inverse  relations 
clearly  suggests  numerical  treatment.  Here,  ir  >  0,  and  K3  is  chosen  so 
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that  it  is  always  real.  From  Eq.  (12A),  in  terms  of 
coordinates, 


(uhjHpwhjh^)^  +(puh2h3)^  (whj) 

+  -  ijtpu^h,)  (uh3)(pwh1h2)^  -  (puh2h3)^  (whj)  ^ 


±  a'e, 


Kj|j  -  a'-2{puh2h3)2 


1  + 


(uhjHpwhjh^  -  (puh2h3  )|3(whi, 

[uh3)(pwh1h2)^  -  (wh1)(puh2h3)^ 


1/2 


(66) 


and 


*■  = 


puh2h3 


e2  (u^Kpwh^)^  +  tpuh^^wh^  _g] 


1  -  %\  (uh3)(pwh1h2)e2  -  (puh^^^wh^  c2  +  e2 


(puh2h3) 


(uh  Kpwh.h 


12^ 


+  (puh,h^)t  (whj) 


1 


1 


\-%\~  (uh3>(pwhlh2)e2  -  ^“Wg^l7  c2  +  g2 


-  p--2(c2  +  e3)(i  -  ifp 


1/2 


(67) 
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which  is  of  the  indicated  form  of  Eqs.  (35)  and  (36)  (the  A^,  A^,  A^ 
functions  have  £3'  u>  v>  w  as  arguments),  therefore,  providing  a 

means  of  construction  of  the  <P  surfaces. 


IV .  THE  TRANSFORMED  SPACE  AND  THE 
RESULTING  SYSTEM  OF  EQUATIONS 


In  the  prior  section,  the  inverse  relations  necessary  for  Eqs.  (5)  through 
(9)  to  be  written  entirely  in  terms  of  (tr,  t,  tt)  coordinates  involved  an 
approximate  inversion  ofEqs.  (62)  through  (65)  'with  u,  v,  w,  p,  held  fixed. 
The  transformed  system  is,  therefore,  determined  to  within  as  close  an 
approximation  as  desired  in  a  neighborhood  of  a  point  in  (cr,  t,  tt)  coordi¬ 
nates  mapped  from  (£^,£->£3)  coordinates  byEq.  (12)  or  (12a).  For  a 
numerical  marching  procedure  involving  finite  difference  methods,  on  the 
basis  of  the  transformed  equations,  the  additional  calculational  error  due 
to  systematic  errors  for  the  approximately  determined  functions  in  small 
regions  may  be  made  to  possess  a  small  effect  in  total  error  by  suitable 
adjustment  of  mesh  size,  etc.  ,  in  the  finite  difference  approach. 

*V»  >X>  r\,  r\,  .  - 

At  a  point  in  R^  at  which  velocities  u,  v,  w  exist,  then  it  follows  that,  with 
hl  =  h1(<r- T>  tt;  u,  v,  w,  p)  ,  h2  =  h2(cr,  t,  it;  u,  v,  w,  p)  ,  h3  =  h3(cr,  t,  it;  u,  v,  w,  p)  , 


hjU 

Mi 


T  -  a 


(uh3)(pwhIh2)^  +  (puh2h3)^  (whj) 
1  *"  * _ 1_  1  1 


pU  2  3  (uh3)(pwhjh2)^  -  (puh2h3)^  (whj) 


'-It  ,  ,  £  '-1  ,  ,  (uh,)(pwh,h,)t  -  (puh,h,),  (wh  ) 

a  h1vpuh2h3  h^a  puh2h3  3  ,r  1  2'?3  r  2  3  1 


Ml 


Ml 


(uh3)(pwh1h2)e 


(whj)(puh2h3)^ 


(68) 
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V 


£ 

V 

-^1 

hl 

Cz+i\ 

+ 


-p  *puh2h3 

(uh3)(pwhjh2)^  +  (puh2h3)^  (whj) 

[M)M 

|jl/2  (uh3)(pwhih2)^2  -  (Puh2h3>e2(wh1) 

P 

~  i  i 

h2wp  -1 

imm: 

iji/2i 

1  h3  [(c2  +  £?)  (i  - «; 

or 

(uh3)(pwh1h2)^  +  (puh^Jg  (whj) 

PVh!  3  puh2h3  (uh3)(pwhjh2)^  -  (pu^h^  (whj) 


(69) 


w 


'v, 


h3UTT/  1  ,  1  9t  ,  1  9<r  \  x  h3VTT/  J_  J_  9T 
hl  V^l  T  8gl  ^  ®V  h2  \  e2  T 

h3WTT  .  !  i  8X  l  0(j-  \ 

h3  '  ^3  T  a^3  ^  a^3' 


+ 


(70) 


where  the  inverse  relations  have  been  used  throughout  in  the  metrics  and 
the  (^,^2>^3 )  system.  The  u,  v,  w  are  then,  respectively,  parallel  to 
the  <r ,  r,  tt  axes . 


If  the  contravariant  form  for  the  velocity  is  assumed,  v  =  Sh  vne  ,  where 

7  n  n  n 

v11  =  v  /h  and  v11  are  the  components  of  a  contravariant  vector,  then  the 
n  n 

transformation  of  components  in  R  and  is  related  by 
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'V. 

where  are  t,  <r,  ir,  respectively,  for  n  =  1,  2,  3;  these  forms  are 
essentially  Eqs.  (68)  through  (70).  In  the  derivation  of  the  covariant 
components  of  Eq.  (45),  in  general, 


f1. 
>  J 


(14b)  or  (72) 


where  f1  are  the  covariant  components,  and  (,  j)  denotes  differentiation  with 
respect  to  In  performing  operations  on  i|>,  0  derivatives  with  Eq.  (14b), 
it  may  be  seen  that  Eq.  (45)  results. 


It  remains  now  to  express  the  original  system  in  terms  of  the  approximate 
system  (for  numerical  analysis  purposes)  in  terms  of  (cr,  t,  it)  coordinates. 
The  form  (Eqs.  5a  through  9a)  is  invariant  in  all  coordinate  systems 
herein.  From  previous  considerations,  then,  the  (cr,  t,  it)  system  in  Rj  is 
of  the  general  form  (excluding  the  energy  equation  for  the  present) 


P  dxJ 


0 


(73),  (74),  (75) 


(pv  )  .  =  0  (76) 

)  1 

where  the  metrics  of  the  transformed  system  are  known  as  combinations 
of  the  metrics  of  the  original  system;  the  velocities  of  the  original 
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system,  written  in  contravariant  form  are  then  simply  contravariantly 
transformable  into  the 


U  V  w 


contravariant  components  in  R 

Here 


In  Eqs .  (73)  through  (77)  the  x1-1  independent  variables  represent  (< r,  t,  it), 
respectively,  for  j  =  1,  2,  3.  In  that  a  finite  difference  approach  has  been 


suggested  throughout,  these  transformed  equations  in  with  known 
boundaries  are  to  be  solved  in  R^;  solutions  in  Rj  are  then  obtained  through 
Eqs,  (62)  through  (70),  and  (12a).  In  Eq.  (73),  p,  p  are  taken  as  scalar 
quantities  and  transform  as  scalars.  The  new  form  of  the  energy 
equation  will  now  be  discussed. 


In  Ref.  1,  f(Ty2)  =  PR"Y 


arises  from  the  integral 


dT 

(-r2)/y 


dy 


(77) 
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On  a  streamline  Eqs.  68  through  70  become 


■  Th, 


v 

h. 


stm. 


stm. 

hl^l 

stm. 

u4l 

ve2 

I  (c2  +  ^)hi  hz(l  -  *2) 


stm. 


'V. 

W 

+  v  1  W 

'■U 

h3 

—  —  II 

stm. 

.(c2+5?)h,  ^ 

stm. 


(78) 


(79) 


(80) 


and  the  u,  v,  w  components  in  may  be  expressed  in  terms  of  deriva¬ 
tives  of  the  independent  variables  by  Eq.  (2  7)  along  a  streamline.  In 
general,  <r  =  ir^,  £2>  £3;  u,  v,  w,  p),  t  =  t(^,  £2>  £3;  u,  v,  w,  p),  and 
tt  =  it(£  ^ .  ^2’  ^3»  u>  v>  w»  P)  and  inverses  were  numerically  determined  in  the 
neighborhood  of  a  point  with  u,  v,  w,  p  fixed  as  £  =  £  j(cr,  t,  tt;  u,  v,  w,  p)  , 

l2  =  ^2(o".  t,  it;  u,  v,  w,  p)  ,  £3  =  £3(<r,  t,  it;  u,  v,  w,  p)  or  alternatively,  as 
=  £  1  (<r ,  t,  ir;u,v,w>p)  ,  £2  =  £2(cr ,  t,  tt;  u,  v,  w,  p) ,  £3  =  £3(o-,t,  n;u,  v,  w,  p) . 
On  a  streamline,  however,  ratios  of  velocities  in  a  particular  coordinate 
frame  are  known  as  ratios  of  derivatives,  as  indicated  for  the  < 
surfaces  by  Eqs.  (40)  through  (43).  Therefore,  with  Eqs.  (28)  and  (29),  the 
two  integrals  representing  stream  surfaces  arise  from 


t  .  -dT/  t 

r  dV^i 


stm. 


-dir/  it 


^ld^l 

c2  + 


stm. 


(81) 
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(82) 


I 


dcr/cr 


2‘  -^1  dei  jl  d^2 

T  - 


c2  +  i\ 


1  -  e: 


-dir/ it 


stm. 


^1  dei 

c2  +  e2 


d£0  d£ , 


stm. 


2  2  2  1 2 
where  >  C  ,  and  <  1  outside  of  the  «  |  stagnation  streamtube 


f\j  Oj 


approximation  surface.  Provided  h1  ,  h^,  and  the  remaining  geometri¬ 
cal  quantities  of  Eqs.  (78)  through  (80)  are  prescribed,  then  velocity  ratios 
on  streamlines  are  determined.  On  a  streamline,  from  Eq,  (12a),  or 
Eqs.  (81)  through  (82), 


const. 


<r  = 


const. 


7\11/2 

(c2  +  6?)(1  ■ 

(83),  (84) 


const. 


^2^3 


IMM] 


1/2 


(85) 


Next,  for  a  particular  tr ,  t,  ir  (starting  from  a  location  on  the  shock  front), 
it  is  then  possible  to  determine  relations  of  the  form  =  ^(cr,  t,  it)  , 

'X. 

^2  =  ^2^°"’  T’ ^  ’  ^3  =  °n  stream^nes  R-  Provided  the  tangent 

streamline  is  next  constructed  over  a  small  distance  from  the  original 
surface,  new  (cr,  t,  it)  values  are  obtained  in  such  a  manner  that  the 
streamlines  in  (<r,  t,  tt)  space,  R,  are  a  priori  determinable.  Then  in 
form,  locally, 

pp’^  =  f[  I^tr,  t,  it),  I^tr  ,  t,  it)]  (86) 
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Although  this  procedure  places  heavy  emphasis  on  numerical  routines  to 

[  r\j  t\,  r\,  r\j  r\^  f\j 

and  surfaces  u/u,  v/v,  w/w,  for  cr ,  t,  it,  h.  ,  h_  ,  h_ 

1  c  5 

specified  fromEqs.  (83)  through  (85),  and  (78)  through  (80)  such  that 
^i  =  ^i^’  T’  1T>U>  v>  w)  =  ^(tr>  T>  v,  w),  etc.  ]  ,  essentially,  these  are 
necessary  steps  for  the  particular  transformation  chosen  (Eq.  12). 


Also,  as  -  constants,  in  Eqs.  (73)  through  (75),  and  for  velocities 

U,v,W’x\i  j/  =  Provided,  next,  the  Eq.  (75)  is  replaced  by  Bernoulli's 
equation,  then  together  with  Eq .  (86),  after  elimination  of  p,  p  terms, 
three  nonlinear  coupled  equations  in  terms  of  u,  v,  w  result,  similar  to 
the  end  result  of  Ref.  1.  These  equations  are  then  subject  to  numerical 
methods  for  solution  in  a  manner  analogous  to  that  of  Ref.  1 . 
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V.  SOME  REMARKS  CONCERNING  NUMERICAL  ANALYSIS  SCHEMES 
FOR  IMPROPERLY  SET  ELLIPTIC  SYSTEMS  OF  EQUATIONS 


The  system  (Eqs.  73  through  76)  and  the  energy  equation  (Eq.  86) 
have  been  implicitly  given;  the  explicit  form  may  be  given  either  in  terms 
of  u,  v,  w  or  u,  v,  w  velocity  components.  In  terms  of  the  latter 
system,  and  u,  v,  w  velocities,  the  system  is  of  the  form 


22 


0  (87) 


where  x1  denotes  <r,  t,  it  for  i  =  1,  2,  3,  respectively. 


In  that  the  sonic  surface  may  be  a  priori  determined  in  the  R^  space, 

'\y 

the  transformation  into  the  Rj  space  follows.  The  difference  scheme, 
starting  from  the  shock,  is  then  required  to  satisfy  conditions  on  the 
given  sonic  line.  In  Eq.  (87),  the  elements  of  the  matrices  are, 
respectively, 


P  = 


<»  + 


~2  I/- 
w  ) 


LL_ 1 

2yf 


1/v-l 


(88) 
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where 


„2  ,-^2  ,  ~2  ,  ~2/|/y  -  1  1  t 

Coo  -  (u  +  v  +  w  u  f 


(89) 


„2  ,  ,  2  , 

^  =  1  + - *-  ,  and 

oo  ,  , ,,  ,2 


( Y  -  DM 


oo 


„1  =  J_ 

"H  y 
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2  ,^2  ,  -^2  ,  ~2x"j  Ty  -  I 
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00 
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oo 
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^2  ~2  ~2."|-  v/v- 1 

(u  +  v  +  w  )  Y  Y 

<~2  ~2  ~2."1- v/v- 1 

(u  +  V  +  w  )  Y  Y 

,~2  ,  ^2  ^2 1- v/v- 1 
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1/7-1 


-  2 
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00 
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-v2 
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,~2  ~2  ^2V 
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The  system  Eq.  (87)  is  then  subject  to  finite  difference  solution  in  the 
manner  of  Ref.  1;  solutions  at  each  point  of  the  mesh  may  be  trans¬ 
formed  into  the  physical  space  by  the  inverse  relations  between 
Rj  and  R^  previously  developed.  The  above  problem  is  then  reduced 
to  simultaneous  solution  of  algebraic  equations  at  points  of  the  mesh. 


To  conclude,  finite  difference  methods  for  problems  of  this  nature  have 
been  discussed  in  Refs.  4  and  5  for  the  plane  case  of  the  Laplace 
equation.  A  discussion  of  Ref.  4  has  been  given  in  Ref.  1;  a  brief 
discussion  of  Ref.  3  will  be  given  here,  in  that  methods  of  this  type  are 
suggestive  of  methods  to  be  used  for  the  problem  in  question. 
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1 2 

For  the  Laplace  equation  u  +  u  =  0  ,  u(0,  y)  =  u(ir,  y)  =  0  , 

_  oo  , 

u(x,  0)  =  <p(x)  ,  u  (x,  0)  =  ip(x)  ,  then  u(x,  y)  n/tt / 2  =  a,  e  y  sin  kx 


oo 


+  X,  bt 


rkv 


k=l 


k=  1 


sin  kx;  in  Cudov's  subsequent  analysis  of  the  Cauchy 


problem,  the  second  term  of  the  solution  is  dropped  throughout  to 


simplify  the  results.  Then  provided  in  0  <  x  <  it  ,  0<y<Y,  with 
rv 


f 


2  2 

u  (x,  y)  dx  <  M  ,  the  norm  on  the  initial  data  is  taken  as 


^  2  2k  Y 

X  ake 

k=l 


1/2 


(90) 


and  the  norm  on  u  is  taken  as 


I  u  I  I  =  sup  \  f  U2(x,  y) 
0<y<YLA 


1  1/2 


dx 


(91) 


then  |  juj  |L  |  jpj  |  q.  In' terms  of  these  norms,  the  effect  of  round-off 
errors  in  various  difference  schemes  are  determined  by  assigning 
approximate  values  to  the  function  u(x,  0)  =  p(x)  on  x  =  x^  =  nh,  with 
h  =  tt/N,  n  =  1,  ■  •  •  ,  N  -  1  and  then  ascertaining  errors  between 
exact  and  approximate  solutions  on  x^  lines  for  0  <  y  <  Y.  The  "size" 


12 

Unique  solutions  that  depend  continuously  on  the  initial  data  under 
certain  restrictions  on  the  class  of  solution  functions  are  considered 
in  Sudov's  paper. 
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of  u(y)  on  the  mesh  is  then  taken  as 


i  I  u(y) 


(92) 


Next  for  u(x,  0)  =  0(x)  discretized  as 
n  =  1,  .  .  .  ,  N  -  1,  then  u(x,  y)  =  n/ Yj 


u(xn>  0)  =  u(xn,  0)  -  6n 

tt  a^eky  sin  kx,  [where 

k=  1 


are 


determined  from  the  u(xn,  0)  boundary  condition]  so  that  the  estimate 
on  x^  lines  becomes  from  (Eq.  92), 


|u(y)  -  u(y)  | 


,<  \^2  M 
h 


e~{N+|)(Y-y)  c  (N.1W 

-  e-4N<Y-y) 


(93) 


where  M  is  a  constant  and 


6  = 


'  N-l 
h  £ 


Now  for  a  system  of 


L  n=  1  J 

Cauchy-Riemann  equations  on  the  uniform  rectangular  net  with  steps 
Ax  =  h,  Ay  =  t  —  T(h)  with  T(h)  -*  0  as  h  -*•  0,  a  norm  on  the  initial  data 
similar  to  (Eq.  90)  is  chosen 


Oh 


rN-i 

EA2  2k Y 
a  e 

Lk=  1 


1/2 


(94) 


where  the  <1,  are  coefficients  formed  from  the  expansion  of  cp  =  </>(nh) 

K  H 

in  terms  of  the  orthogonal  system!  ^Z/tt  sin  k.nhl  ,  k  =  1,  .  .  .  ,  N  -  1. 
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Now  the  difference  scheme  in  t,  h  is  stable  for 


My) I  lMKl  Ml 


0  <  y  <  Y 


0  <  h  <  hQ 


If  the  difference  scheme  has  solutions  of  the  form 
ufmh,  nr)  =  =  sn(k;h)  sin  kmh,  then  the  necessary  condition  of 

stability,  from  (Eq.  95),  is  |  sn(k;h)  |  <  Ke^^,  (0  <  hT  <Y),  with  a 
strong  form,  |  sn(k;n)  |  <  Ke^,  0  <  nT  <  y.  If  next  1  / t  [log  s(k;h)]-*-k 
as  h  -*■  0,  for  0  <  k  <  k^,  then  converges  to  the  exact  solution  in  the 

V 

norm  Eq.  (92).  Cudov  then  proceeds  to  give  seven  difference  schemes 
with  conditions  on  r,  h  such  that  the  strong  stability  condition  is  met, 
and  also  gives  the-error  of  approximation  for  the  Cauchy- Riemann 
systems,  i.  e.  , 


error  of  approx.  =  e  =  O(t)  +  0(h  ) 
stable  for  r/h  =  C  =  arb.  const.  , 
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2  2 

error  of  approx.  =  e  =  0(t  )  +  0(h  ) 
stable  for  -r/h  =  C  =  arb.  const. 


etc.  Now  for  the  error  of  approximation  comparable  to  the  round-off 
error  (from  the  stability  of  solutions  with  respect  to  the  initial  data 
strong  stability  criteria),  the  step  size  h  is  determinable  so  that 
maximum  computational  accuracy  may  be  obtained.  Analysis  of  this 
type  is  necessary  for  the  problem  considered  in  (tr,  t,  it)  space; 
however,  extension  to  nonlinear  problems  presents  considerable 
additional  difficulty  so  that  for  purposes\^of  computation,  the  above 
analysis  serves  only  as  a  model  for  estimation  of  certain  norms  such 
that  the  computational  solution  is  manageable. 
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Transformed 


Coordinate  System 


Fig.  2. 
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